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We take the golden-rule instanton method derived in the previous paper [arXiv:1509.04919] and reformulate 
it using a ring-polymer approach. This gives equations -which can be used to compute the rates of electron- 
transfer reactions in the nonadiabatic (golden-rule) limit numerically within a semiclassical approximation. 
The multidimensional ring-polymer instanton trajectories are obtained efficiently by minimization of the 
action. In this form, comparison with Wolynes’ quantum instanton method [P. G. Wolynes, J. Chem. Phys. 
87, 6559 (1987)] is possible and we show that our semiclassical approach is the steepest-descent limit of this 
method. We discuss advantages and disadvantages of both methods and give examples of where the new 
approach is more accurate. 


I. INTRODUCTION 

In the previous paper, henceforth referred to as Paper 
I,^ we outlined a derivation of a golden-rule instanton 
theory for computing electron-transfer rates in the nona¬ 
diabatic limit. This was based on a time-independent 
methodology using Fermi’s golden rule, which is correct 
in the limit that the electronic coupling is weak. In 
these equations, we substituted the semiclassical limit 
of the Green’s functions describing nuclear dynamics on 
one of two potential-energy surfaces at a given energy. 
A number of steepest-descent integrations led to a for¬ 
mula which defines the rate in terms of the action of an 
imaginary-time periodic orbit, known as the instanton. 

In this paper, we show how this approximate formula¬ 
tion of the rate can be evaluated numerically to treat elec¬ 
tron transfer in large complex systems. We describe how 
the golden-rule instanton trajectory can be discretized, 
allowing it to be located efficiently using multidimen¬ 
sional optimization techniques. This is done using a 
ring-polymer instanton approach similar to that used by 
related methods employing a single Born-Oppenheimer 
surface, including the adiabatic rate^“® as well as tun¬ 
nelling splitting calculations. 

In contrast, early applications of instanton approaches 
employed a method known as “shooting” to locate the 
required instanton trajectory. This method ran classical 
dynamics on the inverted potential-energy surface and 
attempted to choose the correct initial conditions such 
that the trajectory closed into a periodic orbit.Because 
the trajectories are unstable, this approach is inefficient 
and in general limited to treating systems of very few 
dimensions. 

Many alternative methods exist for computing nona¬ 
diabatic rate constants based on a time-dependent 
formulation. These include exact wave function 

calculations^®d7 real-time path-integral calcula¬ 

tions for system-bath models.For more general 
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systems approximate trajectory-based methods have 
been developed^^^^^ including extensions of ring-polymer 
molecular dynamics. 

There are some difficulties with time-dependent meth¬ 
ods however, as the flux correlation functions^^ can be¬ 
come very oscillatory when describing electron transfer. 
Some work towards avoiding these problems has been 
achieved by modifying the correlation function formal¬ 
ism to remove the oscillations, although without affect¬ 
ing the long-time limit which defines the exact rate.^^ 
This simplification was achieved in part by considering 
a time-independent picture, as we have also done in the 
derivation of the golden-rule instanton. 

Although the derivation is very different, we also show 
how our result can be related to Wolynes’ quantum in¬ 
stanton method. This approach uses an approximation 
based on the short-time behaviour of the flux correlation 
function in the nonadiabatic (golden-rule) limit and is 
evaluated using path-integral Monte Carlo. The method 
has been applied to study electron transfers in chemically 
and biologically relevant systems.Our new deriva¬ 
tion of a golden-rule rate offers more insight into the ap¬ 
proximations made by such methods and is in some cases 
more accurate. 

An outline of the paper is as follows. The main results 
from Paper I are summarized in Sec. II, and we show 
how the action integral is discretized and its derivatives 
obtained in Sec. III. We thus obtain a ring-polymer in¬ 
stanton formulation for the electron-transfer rate, which 
is related to Wolynes’ quantum instanton approach in 
Sec. IV. Suggestions for how the instanton approach 
could be applied numerically to complex systems are pre¬ 
sented in Sec. V, which introduces an efficient algorithm 
for locating the instanton trajectories. This is applied to 
an example system in Sec. VI to analyse its convergence 
properties, and Sec. VII concludes the article. 
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II. SUMMARY OF THE GOLDEN-RULE INSTANTON 
APPROACH 


It was shown in Paper I that the instanton relevant to 
the electron-transfer problem is an imaginary-time peri¬ 
odic orbit. This is formed of two trajectories which travel 
on the upside-down reactant, Vb(x), or product, Pl(x), /- 
dimensional potential-energy surfaces. They each bounce 
once and join smoothly together at a point, x^, found on 
the crossing seam, defined by Vb(x) = Vi(x). 

In this paper, we deal only with imaginary-time tra¬ 
jectories and thus depart from the notation of Paper I 
by dropping the bar over imaginary properties. The Eu¬ 
clidean action along one trajectory, either on the reactant 
(n = 0) or product (n = 1) surface, 


Sn = Sn{y-',X-",Tn) = 


r 

1 

9x(r) 

L 

2m 

dr 


+ K(x(t)) 


dr, 


( 1 ) 


where the trajectory, x(r), travels through the classically 
forbidden region from x(0) = x' to x(t„) = x", or equiv¬ 
alently in the opposite direction. A complete periodic 
orbit which runs in imaginary time (3% has the action 

S'(x',x",r) = So{x',x",l3Ti - t) + S'i(x",x',t), (2) 

where r S [0,j3h]. The particular periodic orbit required 
is that which is stationary in x', x" and r. In the following 
all terms are evaluated at this stationary point, at which 
x' =x" =xl. 

The golden-rule instanton method derived in Paper I 
gives a semiclassical approximation to the rate in terms 
of the actions along these trajectories. Two equivalent 
formulae are 
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(4) 


in a few special systems for which the bounce trajecto¬ 
ries and corresponding action is known analytically. In 
order to apply the method to more general problems with 
anharmonic potentials, we will require numerical meth¬ 
ods which are able to locate the instanton trajectory and 
evaluate the action and its derivatives. This is the topic 
addressed in this paper. 


III. DISCRETIZATION SCHEME 


In this section, we show how the action integral, 
Eq. (1), can be defined from a discretized form of an 
imaginary-time trajectory. This is based on the ring- 
polymer instanton method,^ which has been successfully 
used in adiabatic, single-surface, rate calculations®’® as 
well as the evaluation of tunnelling splittings.It re¬ 
lies on the fact that a classical trajectory is known to 
give a stationary value of the action, with respect to any 
deviation along its length except at the end points.®®’^® 

We also describe how second derivatives of the action 
can be evaluated directly without resorting to taking fi¬ 
nite differences between instantons optimized under var¬ 
ious conditions. The approach we use for this follows 
closely the method of implicit differentiation described 
in Ref. 8, which we extend to obtain all the derivatives 
required for the golden-rule instanton method. 

According to our golden-rule approach,® we only need 
to study the dynamics on one of the two potential-energy 
surfaces at any time. This section would thus also be 
directly applicable to single-surface reactions, simply by 
dropping the subscript n. 

We consider an imaginary-time pathway of length r„ 
between the points x' = Xg and x" = X 7 v„, which passes 
through the intermediate points {xi,... ,XAr,^_i} at a set 
of discrete times. The imaginary-time intervals between 
each point are = CiTn, with j G {1,..., A^„} such that 
each Ci G [0,1] and = 4- T4ie velocity along a 

given pathway at these times is given by jx^ — Xi_i|/eiT„ 
and the action by 


where the van-Vleck prefactor for a trajectory is given by 


Cn 


dx'dx" 


(5) 


and the other prefactors are 
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These formulae for the golden-rule instanton method 
were used in Paper I to obtain the rate of electron transfer 


where the first term originates from a trapezium-rule in¬ 
tegration of the kinetic energy along the pathway, and 
the second of the potential energy. This is the general 
form allowing for uneven imaginary time intervals®" and 
would simplify to the usual case with = l/Nn- Note 
that here an open-ended pathway is described such that 
no cyclic indices are implied. 

The points x^ for i G A^„ — 1}, which give the 

coordinates along the classical trajectory, are those which 
































3 


give a stationary value of Sn, i.e. those which solve 

(Ci + ei+l)T^ dVr 


x» - Xi_i X, - Xj+I 

e* ei+i 


2m 


i9x, 


= 0. (9) = 


differentiation by gives 
to(x'-Xi) eiT„ 


+ 


The action along the trajectory is therefore 
5 „(x',x",t„) = hmAr^_>oo *5„(x',x",r„), where 

5„(x',x",t„) = S'„(xo,Xi,... ,XAr„;r„) and Xq = x', 

XAf„ = x". 

In fact the dominant classical trajectory between 
two end points in a given imaginary time will be the 
global minimum of Eq. (8) with respect to the inter¬ 
mediate points. This can be obtained by employing 
a multidimensional optimization routine such as the 
limited memory Broyden-Fletcher-Goldfarb-Shanno (1- 
BFGS) algorithm^^ in the same way as is done for 
tunnelling splitting calculations.However, the end 
points, x^, required for the instanton method are not 
in general known a priori, so the instanton trajectories 
cannot be obtained in this way. We discuss optimiza¬ 
tion methods which do not require knowledge of the end 
points in Sec. IV. 

Differentiating Eq. (9) by the end points, x' or x", gives 
equations which can be written in the following form: 
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where the elements of the doubly indexed matrix J are 
dehned by 
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Differentiating again, we obtain the second derivatives: 
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Again the indices are not cyclic, i.e. the matrix is banded 
with bandwidth /. This definition is equivalent to that 
given in Ref. 8 when the time-steps are equal. Equations 
(10) can be solved numerically for the derivatives of x 
using standard linear-algebra routines. Note that these 
partial derivatives imply that t„ and one end point are 
kept fixed while the rest of the pathway is allowed to 
re-optimize itself as the other end point varies. 

Other terms are found by differentiating Eq. (9) by t„ 
to give 


Partial derivatives of Sn are approximated by these for¬ 
mulae, which become exact in the A4 —>■ oo limit. As¬ 
suming that the instanton trajectories have already been 
found, these derivatives can be applied in the prefactor of 
Eq. (3), using Eq. (7) and Eq. (2), to give the golden-rule 
instanton rate. 

In contrast to standard approaches where the eigenval¬ 
ues of a iV/ X Nf matrix are required for the prefactor, 
the most difficult task in this approach is the solution of 
the linear equations, Eqs. (10) and (12). Because J is the 
Hessian matrix about the minimum pathway, it is posi¬ 
tive definite, and the equations can be solved efficiently 
using a Gholesky decomposition, taking advantage of the 
banded nature of the matrix. 


Nii-l f 

E/ E/ 

i'=l i'=l 


dxi' jf 
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Nnm 9ij^- ’ 


( 12 ) 


which we solve for . 

Using the fact that Sn is stationary with respect to 


This approach is not limited to the current application 
but may also significantly improve the efficiency of other 
instanton methods, for which the diagonalization can be 
a considerably time-consuming task for high-dimensional 
systems. We shall discuss the use of such an approach to 
improve the efficiency of adiabatic rate calculations in a 
forthcoming paper. 
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FIG. 1. Schematic showing the ring-polymer beads which 
discretize the instanton for the case of No = 6 and A^i = 4. 
Those on the left coloured in blue have the electronic configu¬ 
ration of the reactant state |0) and those on the right in red of 
the product |1). Beads No and N are located at the hopping 
point X* and contribute to the action of both trajectories. 


IV. RING-POLYMER INSTANTON FORMULATION 


So far, we have only dealt with open-ended trajectories, 
whose end points are as yet unknown. In this section, we 
extend this methodology to obtain the pathway for the 
total periodic orbit. This orbit is simply the combination 
of a trajectory on the reactant surface with another on 
the product surface and has total imaginary time f3h. 

We divide up the total orbit into N segments, with the 
first Nq on |0) and the remaining 7Vi = N — Nq on |1) 
as in Fig. 1. Equal time-step intervals, = 1/Nn, will 
be assumed here but other choices may slightly improve 
efficiency. ' 

There is a special case that the time intervals on 
both trajectories are equal, all with length /3jv7i, where 
= P/N. This can only be obtained in practice if the 
imaginary times along each trajectory, are 

known a priori. Such cases arise for example if the reac¬ 
tion is symmetric, where the stationary value is known 
to be T = fiTi/2, or if the instanton has already been 
obtained by an alternative method, such as those intro¬ 
duced in Sec. V. 

In this case, we have a formulation similar to path- 
integral^^ and ring-polymer molecular dynamics, 
which were obtained from a discretization of the quan¬ 
tum Boltzmann operator. The resulting set of N co¬ 
ordinates are called beads and, via a quantum-classical 
correspondence,^® are equivalent to a ring polymer 
of classical particles connected together by harmonic 
springs. 

It is a good idea to use the Wbead steepest-descent 
approximation to the reactant partition function,^® 


phujj 


/ r n^~ 1-1 

Zq = 2 sinh ' 

i=i 

2 . , P]\f?iujj 




sinh 


(15) 

(16) 


as this is known to benefit from a cancellation of errors 
with the iV-bead instanton calculation and improve con¬ 
vergence of the rate.® Here ujj are the normal-mode fre¬ 
quencies at the minimum of Vo(><); if there are translation 


or rotational modes, the formula should be modified ap¬ 
propriately. 

The total action along the two joined pathways is given 
by 


I3nUn{x) = S'o(xAr,Xi, . . . ,XNg-, Nq/Sn^) 

+ <S'i(xAro,. ■. ,X7v; IVi/SatTi), (17) 
such that the iV-bead ring-polymer potential is 

Un{x) = 

i=i ^Pn'^ 

No-l 

+ ^Vb(xiv) + ^ Vb(x2) + |VJ)(xAro) 

N-1 

+ I^i(XA^o) + ( 18 ) 

i=No + l 

The positions of each bead are given by a; = {xi, ..., xjv}, 
and cyclic indices are implied such that Xq = xjy. This 
function can be minimized with respect to the positions 
of all beads to obtain the coordinates x = {xi,... ,xn} 
along both trajectories simultaneously. The hopping 
point is then identified as x^ = xjvq = xjv and the ac¬ 
tion a.s S = PnUn, where Un = Um{x). However, r, or 
equivalently the ratio Nq/Ni, is yet to be specified. It 
will therefore be necessary to compute the instanton for 
numerous values of r to find the stationary point of Un 
with respect to r. 

We now introduce the quantum instanton approach of 
Wolynes.®^ This method was derived using a steepest- 
descent evaluation of the time integral over the exact 
flux-flux correlation function within the golden-rule ap¬ 
proximation and gives 

= (19) 

^ ^-Nf j g-/3«C/„(x) (20) 

where the prefactor is A = 27r/3Ar?i^/m. It is here as¬ 

sumed for simplicity, that the electronic coupling, A, is 
approximately constant, although the formulation could 
be generalized without affecting our findings. In practice 
the integrals are computed using a discrete path-integral 
Monte Carlo simulation, and r = Ai/3jv7i is chosen in 
the range [0,/3?i] such that ^ = 0. Taking the second 
derivative of Eq. (20) gives^^ 

/ ^-(x')bC(x")e-^«^~(=") da:, 

(21) 

where H-(x) = Vo(x) — Vi(x). The second derivative of 
(j) is negative and thus corresponds to a stationary point 
which is a maximum along r. 
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The derivation is similar in spirit to that used to obtain 
the quantum instanton approach for Born-Oppenheimer 
systems described in Ref. 48, as it also employs a 
steepest-descent integration along the real-time coordi¬ 
nate of a flux-flux correlation function. The single¬ 
surface quantum instanton approach is however not a 
semiclassical approximation in the sense that it gives the 
correct leading order of Ti. This is most easily seen from 
the fact that it does not reproduce correct results for 
a free-particle or in the classical limit.Wolynes’ for¬ 
mula, Eq. (19), is also not exact in the high-temperature 
limit in general. However, it is known that it reproduces 
the stationary-phase approximation^'^ for the golden- 
rule rate of a spin-boson system and hence also Marcus 
theory,®'^ which is the correct result for this system in the 
classical limit. 

To show the link between the quantum and semiclas¬ 
sical instanton methods, we perform a steepest-descent 
approximation to Eq. (20) in two steps, reserving the in¬ 
tegrals over beads assigned to x' and x" until after all 
others. This gives 


^-<j>{T)/n 



X 



g-'So(x',x",/3fi-T)/?i-Si(x",x',T)/?i (Jx'dx^^ 


( 22 ) 


C'oC'i -s/n 
C 


(23) 


Using this definition of (/)(t) in Eq. (19) gives a rate which 
is not in general equal to that of classical golden-rule 
transition-state theory.'^^ This is most easily seen from 
the example of the transfer from a harmonic oscillator 
to an anharmonic product state, such as the system 
discussed in Ref. 33. As shown in Paper I, the high- 
temperature golden-rule instanton rate gives the exact 
classical golden-rule transition-state theory limit for this 
one-dimensional system, 

fcci.TST^o = ^\l~^ J ^[Vb(a;) - Vi{x)\ dx, 

(26) 

whereas Eq. (25) noticeably does not include a delta func¬ 
tion constraining the integral to the crossing seam and 
thus gives an incorrect result. 

This is at first sight surprising, as one would naively 
assume that the steepest-descent approximation reduces 
the accuracy of the result. The reason for the discrepancy 
is that the two methods are based on different approxi¬ 
mations. This example makes it clear that, at least for 
certain problems, a more accurate quantum rate theory 
is obtained from semiclassical considerations than from 
Gaussian approximations to the flux correlation function. 

The link between the semiclassical and quantum in¬ 
stanton approaches also suggests that another method 
could be used to compute the golden-rule instanton rate, 
where the steepest-descent integration is taken over all 
ring-polymer beads simultaneously giving 


where Jn is dehned in Eq. (11) with Ci = 1/Nn and we 
have used the following result from Ref. 8: 
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^—0nUn 


(27) 


Also, within the steepest-descent approximation, « 

12 o 

and therefore, taking these semiclassical limits in 
Wolynes’ formula, Eq. (19), reproduces the golden-rule 
instanton rate, Eq. (4). This shows a strong link between 
the semiclassical instanton theory presented in this pa¬ 
per and the quantum instanton approach—the former 
is a steepest-descent approximation to the latter. The 
quantum instanton approach has a great advantage over 
the semiclassical instanton method, which is that it can 
also treat liquid systems, where many minima exist on 
the ring-polymer potential surface. 

Note however that Wolynes’ quantum instanton is not 
always more accurate than the semiclassical instanton. 
In the high-temperature limit, the ring-polymer beads 
collapse and Eq. (20) reduces to give an integral over the 
centroid mode. 


lim 

/3^o 



g-[(/3?i-r)Vo(x)+-ryi(x)]/?i 

( 25 ) 


where V^Uat is the Hessian matrix found by differentiat¬ 
ing Eq. (18) by all bead positions x^ and is evaluated at 
the instanton geometry, x. 

Because the steepest-descent integrals are evaluated at 
the hopping point where U_(xt) = 0, we have to consider 
a higher-order term for our semiclassical approximation 
of Eq. (21). This is 


h dr^ 


A-Nf r ptA(xt) 

~ir isD 



dV.jxi) 

dxi 



g-/3„C/«(cc) 


(28) 


We evaluate the integral using a second-order expansion 
of Un{x) about the ring-polymer instanton orbit, which 
gives. 


dr^ 


1 au-(xt) 

h dxi 




dV-jxi) 

dxi 


(29) 


where only the / x / block corresponding to rows for bead 
Nq and columns for bead N is required from the inverse 
of the full Hessian. The golden-rule instanton rate in 
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ring-polymer form is thus 


kscZo = V2TTh, 





V^Un 



1 

2 

^—I^nUn 


(30) 


This formula gives the same result as Eqs. (3) and (4) in 
the N ^ oo limit. 

Note that all eigenvalues of the Hessian are positive. 
This is therefore a more straightforward derivation than 
is achieved using the Imf approach,where the in- 
stanton has a negative eigenvalue, which has its sign re¬ 
versed, and a zero-mode which has to be integrated out 
analytically. 

As in the adiabatic, single-surface, case,^ this ring- 
polymer instanton approach provides a computationally 
tractable way to obtain the reaction rate of a multidi¬ 
mensional system. However, it would be necessary in 
general to optimize Eq. (18) many times to find the value 
of r which gives a maximum value of Un- In Sec. V, we 
shall propose alternative methods which obtain r auto¬ 
matically from a single optimization and may therefore 
be found to be more efficient in practical applications. 


V. NUMERICAL EVALUATION 

In this section we present two methods which we 
suggest could be used to evaluate semiclassical golden- 
rule rates in complex multidimensional systems. It 
may also be possible to implement similar schemes 
for locating other instantons more efficiently, including 
those for adiabatic rate theory^ and tunnelling splitting 
calculations.Applications of the methods to such 
systems will be explored in future work. 

In Sec. IV, we discussed an approach similar to that 
used for adiabatic instantons, where the imaginary time 
of each trajectory is chosen before the ring-polymer in¬ 
stanton is optimized. Here we present two alterna¬ 
tive methods which optimize all unknown variables si¬ 
multaneously. The first is based on a Lagrangian for¬ 
malism with equal time-steps and the second uses the 
Hamilton-Jacobi abbreviated action with evenly spaced 
ring-polymer beads. 

Note that the symmetry of the instanton pathway can 
be used to reduce the number of independent coordinates 
to Af/2 -f 1.^ It is known that the instanton must fol¬ 
low the same pathway in both directions of its periodic 
orbit, such that we only need to optimize two shorter 
open-ended trajectories, each with one end at the hop¬ 
ping point and the other at a turning point. In both 
cases, we will employ the bead ordering given in Fig. 1 
and assume that Nq and Ni are always chosen to be 
even. There is a symmetry equivalence between the top 
and bottom rows such that when the pathway is op¬ 
timized, x^o/ 2 _i = ^No/ 2 +i for i G {l,...,No/2} and 
^No+Ni/ 2 +i = ^No+Ni/ 2 -i fo^ * G {Ij • ■ ■ ;-^l/2}. The 


beads at the turning points, >^No /2 and ^No+Ni/ 2 , are in¬ 
dependent. 


A. Lagrangian formalism 

The Lagrangian formalism defines classical trajecto¬ 
ries according to the elapsed time. It was used to de¬ 
fine the standard ring-polymer instanton approach for 
single-surface systems with equal time-steps.As in 
Sec. IV, we again separate each trajectory into Nn equal 
imaginary-time intervals, i.e. with et = 1/Nn. However 
in contrast to the previous approach, Eq. (17), the reac¬ 
tant trajectory may have a different time step from the 
product trajectory. The total discretized action is 

S{x,t) = 2So{y^No/2, ■ ■ ■ ,^No',\{PkL - t)) 

-I-2S'i(xjVq, . . . ,XjVo+Ari/2; 5'^)) (31) 

where due to the forementioned symmetry we have taken 
twice the action along each pathway from the turn¬ 
ing point to the hopping point in half the imaginary 
time. The classical imaginary-time periodic orbit can 
be found as the first-order saddle point of this func¬ 
tion with respect to the independent bead coordinates 
X = {xjVo/ 2 ) • ■ •) XAfo-i-Afi/ 2 } T simultaneously. The 
other half of the instanton orbit is given by symmetry. 

Saddle-point optimization algorithms have been well 
studied in the pursuit of locating instantons,®’®^ in most 
cases a Hessian-based quasi-Newton method being ap¬ 
propriate. The value of the optimized function gives the 
required total action S in the A^-bead approximation and 
the imaginary times tq = pti — t and ri = r. In the 
N ^ oo limit, this result is in principle independent of 
the choice of the ratio Nq/Ni, although an intelligent sug¬ 
gestion would be tq/Nq « n/A^i to make all time-steps 
approximately equal. 

In this way, it is possible to evaluate Eq. (3) numer¬ 
ically using this ring-polymer instanton approach and 
converge the results obtained with respect to N. We 
identify x' = xjv and x" =XNg, both of whose optimized 
positions tend to x^ in the N ^ oo limit. Derivatives of 
the total action S are given as sums or differences of the 
derivatives of Sq and Si defined in Sec. HI. Note that 
here the full trajectory, from x' to x" is required and not 
just the trajectory to the turning point. 

However, this approach requires a saddle-point opti¬ 
mization which is often more difficult than a minimiza¬ 
tion. In Sec. VB, we describe an alternative method to 
locate the instanton trajectories and evaluate their ac¬ 
tions based on a potentially simpler algorithm. 


B. Hamilton-Jacobi formalism 

A significant feature of the derivation presented in this 
paper is that the energy of the two trajectories must be 
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equal. It would therefore be natural to locate the instan- 
ton trajectory under this constraint rather than directly 
attempting to find the stationary value of the imaginary 
time T. To this end, we will employ a Hamilton-Jacobi 
definition for the action along two discretized pathways 
of Nq and Ni ring-polymer beads with the same energy 
for each trajectory. 

We should take care when computing the discretized 
abbreviated action, as a naive implementation using the 
trapezium rule to approximate Wn would give a function 
with infinite derivatives at the turning points. We there¬ 
fore propose the following functional form to compute the 
abbreviated action along one pathway with energy E: 

W„(xo, . ■ • ,Xn^/2',E) 


Pu{x) 

where between each bead we have used the analytical 
expression for the abbreviated action in a linear poten¬ 
tial, and the factor of two accounts for the return journey 
of the trajectory. The absolute value of the momentum 
is taken such that the function returns real values even 
when beads stray into the classically allowed region. This 
ensures that the function is smooth and well-defined ev¬ 
erywhere as is required by most optimization routines. 
The final optimized pathway should however lie entirely 
in the classically forbidden region. This requirement can 
be easily checked. 

In this formulation, it is necessary to ensure that the 
beads remain evenly spaced without biasing the instan- 
ton pathway. The simplest way to achieve this is to in¬ 
clude a penalty function. 


N„j2 


= Xn E ~ 

i—1 

(35) 

5xi = |xj - Xi_i| 

(36) 

2 

^ E 

i—1 

(37) 


This type of approach has been applied successfully to lo¬ 
cate folding pathways in proteins.®^ However, alternative 
methods based on generalizations of the nudged-elastic- 
band algorithm avoid using penalty functions and may 
be more efficient.The value of the scalar Xn should not 
affect the result of a converged optimization and can be 
chosen by the user to maximize efficiency. 

As in all optimization problems, a good initial guess 
is required to ensure fast convergence to the global min¬ 
imum. Instanton optimizations are best performed in 
stages with increasing numbers of beads and decreasing 
temperatures.^*^ An equally spaced straight line normal 


Nr^/2 


= 2 E 




irnKj 


= v'2m|H(x) - E\ 
kn(Xz) l) 


+ En 

(32) 

(33) 

(34) 


to the crossing seam provides a reasonable starting point 
at high temperatures. 

The imaginary time intervals between each bead are 
evaluated from the derivative of the abbreviated action 
with respect to energy as 


Sn = 


7*n(Xi —l) PniX-i') 

Ki 


(38) 


Thus the total imaginary time along each trajectory is 
Tn = “2 6Ti and we define r = j3h /(1 -|- tq/ti). 

Classical trajectories could be located by optimizing 
the abbreviated action Eq. (32) for a given energy. This 
approach would give the microcanonical instanton rates 
discussed in Paper I. However, it is the thermal rate 
which is of most interest, for which the value of E is not 
known a priori. We therefore use the value of the full 
action in the Hamilton-Jacobi picture. 


S{x,E) = Wo{xn„/ 2 ,...,xn^;E) 

-|-ITl(xArg, . . . ,XjVo+Ari/2; E) -I-/3A. (39) 

This function is minimized with respect to the inde¬ 
pendent beads x = {xno/ 2 i ■ ■ ■ EN 0 +N 1 / 2 } and energy 
simultaneously under the constraint that the pathways 
terminate at a turning point, i.e. Vb(xAro/2) = E and 
^i{'>^No+Ni/ 2 ) = E. Constrained optimization methods 
such as sequential least squares programming are ideal 
for this task. 

This Hamilton-Jacobi approach to locating instan- 
tons has significant advantages over the standard ring- 
polymer instanton approach, where the beads tend to 
accumulate near turning points.^ By forcing the beads to 
be evenly spaced along each trajectory, we expect that 
fewer beads will be required to converge the action inte¬ 
gral. The convergence is further improved by using the 
scheme based on the analytic result for linear potentials. 
Another advantage is that the standard instanton-finding 
methods employ a saddle-point search,^ whereas the new 
approach requires only a minimization. It is usually far 
less computationally demanding to locate the latter type 
of stationary point. 

However, it is known that the evenly spaced pathway 
does not give good estimates for the instanton prefactor, ^ 
even when N is large enough to converge the action 
to a high accuracy. This was also confirmed by our 
own numerical tests, employing the formulae in Sec. HI 
with Eq. (38). It seems that the ring-polymer instanton 
methods described in Secs. IV and V A with equal time- 
steps is better for computing the derivatives whereas this 
Hamilton-Jacobi method with evenly spaced beads is bet¬ 
ter for estimating the action. 

We therefore propose that the following combination 
of the methods presented above is used for computing 
the rate: 

• The Hamilton-Jacobi method can be used to lo¬ 
cate the instanton pathway and find the stationary 
value of r. We also take the action, S, from this 
calculation. 










• Using cubic spline interpolation along the imag¬ 
inary time coordinate,the two trajectories are 
modified to give equal time-steps along each trajec¬ 
tory, and the resulting pathway minimized, keeping 
T fixed. 

• The remaining beads in the two bounce trajecto¬ 
ries are obtained by symmetry and the derivatives 
of the actions, Sq and Si, computed using the for¬ 
mulae in Sec. III. 

• The rate constant can be then be evaluated using 
Eq. (3). 

In Sec. VI, we apply this combined method and com¬ 
pare it with the saddle-point search of Sec. V A to a 
model problem. 


TABLE I. Results for the two numerical methods, Lagrangian 
and Hamilton-Jacobi (Ham-Jac), described in Secs. VA and 
VB. In both cases we take Ni/Nq ~ 0.3, although ensuring 
that No and Ni are even. The semiclassical instanton results 
are given in the final row, computed using formulae from Pa¬ 
per I. 


N 


Lagrangian 

Ham-Jac 

s/h 

r/pn 

ksc/kci 

S/h 

ksc/kci 

8 

6.558 

0.3248 

23.2 

6.152 

31.1 

16 

6.179 

0.3163 

33.4 

6.051 

36.5 

32 

6.058 

0.3131 

36.3 

6.020 

36.0 

64 

6.022 

0.3119 

36.1 

6.013 

36.2 

128 

6.013 

0.3117 

36.2 

6.012 

36.3 

256 

6.012 

0.3116 

36.3 

6.011 

36.3 

00 

6.011 

0.3116 

36.3 

6.011 

36.3 


VI. APPLICATION TO A MODEL SYSTEM 


We consider a numerical application of the golden- 
rule instanton method to a spin-boson model of electron 
transfer.®*^’®® Note that the methods are also directly ap¬ 
plicable to anharmonic systems, but here we intend to 
compare with the exact results, which are easily avail¬ 
able only for integrable systems. 

The spin-boson model was defined in Paper I and we 
use the same notation here with parameters chosen to de¬ 
scribe condensed-phase electron transfer at typical con¬ 
ditions. The temperature is T = 300 K, and the spectral 
density of the bath has Debye form J{oj) = f , with 

the characteristic frequency Wc = 500 cm“^, and reorga¬ 
nization energy A = 40kcal/mol. The spectral density is 
discretized with / = 12 bath modes using®®’^^ 


ojj = Wc tan 


0 - 1 ) ^ 

2 / 



(40) 

(41) 


one-dimensional maximization, as 36.3 kd- This is close 
to the quantum golden-rule rate, which was found to be 
36.6 kci by numerical integration. Here, as was also ob¬ 
served in Ref. 35, nuclear tunnelling has a significant ef¬ 
fect on the rate. 

The two numerical approaches outlined in Sec. V were 
applied to the model for various numbers of ring-polymer 
beads. In each case, the starting point for new instanton 
searches was given by a spline interpolation^^ of the tra¬ 
jectories from previous optimizations with fewer beads. 
The results are given in Table I. 

As expected, the rates obtained by both numerical 
methods tend to the semiclassical results in the large N 
limit. The Hamilton-Jacobi formulation is found to give 
better estimates of S than the Lagrangian formulation 
for the same number of beads. Using the combined ap¬ 
proach in which this action is used alongside the deriva¬ 
tives found from an optimized instanton with equal time- 
steps, requires in each case about half as many beads for 
the same error in the rate. This would lead to a signifi¬ 
cant advantage when treating more complex systems. 


Vll. CONCLUSIONS 


where j S We include a bias to products of 

e = 10 kcal/mol. 

The electronic coupling. A, is constant, but for the 
purposes of generality we do not specify its value. It 
must of course be small enough that the golden-rule ap¬ 
proximation is valid. Results are presented relative to 
the classical rate such that they are dimensionless and 
do not depend on A. It was found that 12 bath modes 
are enough to converge the ratio to less than 2%. 

For this model, the classical rate is given by Marcus 
theory as®^ 


kc\ 


. /^p-/3(A-e)V4A 

ny X 


(42) 


Formulae presented in Paper I give the semiclassical 
golden-rule rate, fcsc, with r obtained numerically by a 


In this paper we have described a ring-polymer formu¬ 
lation of the golden-rule instanton approach derived in 
Paper I.^ This formulation is amenable to efficient nu¬ 
merical evaluation and we have suggested two methods 
for its computation. 

The method based on the Hamilton-Jacobi formalism 
appears to be more efficient at obtaining the instanton 
trajectory and its action. This approach forces the energy 
along both instanton trajectories to be equal, which is a 
fundamental aspect of our time-independent derivation. 
Similar approaches may also prove efficient for locating 
instantons used in other calculations, such as adiabatic 
rate theory and tunnelling splitting calculations. 

The ring-polymer instanton was shown to be equiva¬ 
lent to a steepest-descent evaluation of Wolynes’ quan¬ 
tum instanton approach,thus providing a strong 
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link between the two methods. Quantum instanton 
approaches'*® employ a Gaussian approximation to the 
flux-flux correlation function whose short-time behaviour 
is computed using exact path-integral methods. No¬ 
table deviations from Gaussian behaviour occur even 
for the simplest problem of free-particle propagation^® 
and it seems that the flux-flux correlation function can¬ 
not be assumed to be Gaussian if a rate theory is re¬ 
quired which gives a good approximation to the high- 
temperature limit. The golden-rule instanton method 
does not however suffer from these problems. 

All instanton methods will fail when the potential- 
energy surfaces exhibit oscillations, as occurs with liq¬ 
uids, such that many minima appear on the ring-polymer 
surface. In this case, the steepest-descent integrals em¬ 
ployed in the instanton derivation are not valid and path- 
integral sampling methods such as Wolynes’ approach, 
Eq. (19), are necessary. However, for systems where the 
environment is not fluxional, such as in solids®® or cer¬ 
tain gas-phase molecules, the instanton approach may be 
more accurate as well as much more efflcient. 
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